clear
capture log close#delimit ;set more off;version 11.0;
*************************************************************;
*	Replication files for Conrad & Ritter (2013)	    *;
*************************************************************;
use "ConradRitterAnnualData19June12.dta", clear;
sort cowcode year;
**Generating Interactions**;
gen JEJS=linzerstatonJI*jobsecurity_cheibub;

capture program drop biselect;
program define biselect;
*version 7;
    args lnf theta1 theta2 theta3 zrhoS zrhoN;

	tempvar rhoN rhoS;
	qui gen double `rhoN' = (exp(2*`zrhoN')-1) / (exp(2*`zrhoN')+1); 
	qui gen double `rhoS' = (exp(2*`zrhoS')-1) / (exp(2*`zrhoS')+1);

*Selection parameter (rho) from eq2=rhoS, from eq3=rhoN;
      quietly replace `lnf' = ln(($ML_y1==1)*binorm(`theta1',$ML_y2*`theta2', $ML_y2*`rhoS')+($ML_y1==0)*binorm(-`theta1',$ML_y3*`theta3', -$ML_y3*`rhoN'));
          
end;

set matsize 150;
capture drop y1;
capture drop y2;
capture drop y3;
capture drop y1n;

gen y1=catr;
gen y1n=-y1;
gen y2=2*CIRIhightort -1;
gen y3=y2;
	
local ivar1 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS IO;
local ivar2 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;
local ivar3 BanksMobilization linzerstatonJI jobsecurity_cheibub JEJS;

prob y1 `ivar1';
local ll1=e(ll);
qui matrix i1 = e(b);
qui matrix coleq i1 = eq1:;

prob CIRIhightort `ivar2'; 
local ll2=e(ll);
qui matrix i2 = e(b);
qui matrix coleq i2 = eq2:;

prob CIRIhightort `ivar3' if y1==0;
local ll3=e(ll);
qui matrix i3 = e(b);
qui matrix coleq i3 = eq3:;
 
heckprob CIRIhightort `ivar2', select (y1 = `ivar1'); 
qui matrix i4 = e(b);
qui matrix rS=i4[1,"athrho:_cons"];
qui matrix coleq rS= eq4:;

heckprob CIRIhightort `ivar3', select (y1n = `ivar1');
qui matrix i5 = e(b);
qui matrix rN=-i5[1,"athrho:_cons"];
qui matrix coleq rN= eq5:;

qui matrix initvec = i1,i2,i3,rS,rN;

ml model lf biselect (y1 = `ivar1') (y2 = `ivar2') (y3 = `ivar3') () ();
ml search;
ml maximize, difficult;

est store est1;

disp (exp(2*[eq4]_cons)-1) / (exp(2*[eq4]_cons)+1);
disp (exp(2*[eq5]_cons)-1) / (exp(2*[eq5]_cons)+1);
test [eq2=eq3]:_cons;

*************************************************************;

**Calculating predicted probabilities**;

	preserve;
	set seed 10101;	drawnorm CR_b1-CR_b18, n(10000) means(e(b)) cov(e(V)) clear;
	save simulated_betas, replace;
	restore;
	merge using simulated_betas;
	drop _merge;
	postutil clear;	postfile mypost prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi             using sim , replace;            noisily display "start";            	local a=0 ;	while `a' <= 1.0 { ;    {;
    
	*Non-signatory means/modes;
	tab BanksMobilization if catr==0;
	sum linzerstatonJI if catr==0;
	sum IO if catr==0;
	scalar h_BanksMobilization=0;
	scalar h_ls_ji=0.381;
	scalar h_IOnonsig=3;
	scalar h_Constant=1;
	
gen sel = 
	CR_b1*h_BanksMobilization + CR_b2*h_ls_ji + CR_b3*`a' + CR_b4*h_ls_ji*`a' + CR_b5*h_IOnonsig + CR_b6*h_Constant;
sum sel;

gen psel = normprob(sel);
sum psel;
centile psel, centile(2.5,97.5);

gen outcome1a =
	CR_b12*h_BanksMobilization + CR_b13*h_ls_ji + CR_b14*`a' + CR_b15*h_ls_ji*`a' + 
	CR_b16*h_Constant;

gen prob0 = normprob(outcome1a)*(1-psel);
sum prob0;
centile prob0, centile(2.5,97.5);

gen outcome1b =
	CR_b7*h_BanksMobilization + CR_b8*h_ls_ji + CR_b9*`a' + CR_b10*h_ls_ji*`a' + 
	CR_b11*h_Constant;

gen prob1 = normprob(outcome1b)*(1-psel);
sum prob1;
centile prob1, centile(2.5,97.5);

gen diff=prob1-prob0;        egen probhat0=mean(prob0);    egen probhat1=mean(prob1);    egen diffhat=mean(diff);        tempname prob_hat0 lo0 hi0 prob_hat1 lo1 hi1 diff_hat diff_lo diff_hi ;      _pctile prob0, p(2.5,97.5) ;    scalar `lo0' = r(r1);    scalar `hi0' = r(r2);          _pctile prob1, p(2.5,97.5);    scalar `lo1'= r(r1);    scalar `hi1'= r(r2);          _pctile diff, p(2.5,97.5);    scalar `diff_lo'= r(r1);    scalar `diff_hi'= r(r2);      scalar `prob_hat0'=probhat0;    scalar `prob_hat1'=probhat1;    scalar `diff_hat'=diffhat;        post mypost (`prob_hat0') (`lo0') (`hi0') (`prob_hat1') (`lo1') (`hi1')                 (`diff_hat') (`diff_lo') (`diff_hi') ;    };      
        drop sel outcome1a outcome1b psel prob0 prob1 diff probhat0 probhat1 diffhat ;        local a=`a'+ 0.01 ;        display "." _c;        } ;	postclose mypost;                                  	use sim, clear;	gen MV = (_n-1)/100;	graph twoway line diff_hat MV, clwidth(medium) clcolor(blue) clcolor(black) clpattern(solid)        ||   line diff_lo  MV, clpattern(dash) clwidth(thin) clcolor(black)        ||   line diff_hi  MV, clpattern(dash) clwidth(thin) clcolor(black)
        ||   ,               xlabel(, labsize(4))             ylabel(, labsize(4))            yscale(noline)            xscale(noline)            yline(0)            legend(off)            title("", size(3))            subtitle("" "" " ", size(3))            xtitle("Executive Job Security", size(3)   )            ytitle("", size(3))            xsca(titlegap(2))            ysca(titlegap(2))            scheme(s2mono) graphregion(fcolor(white));
*************************************************************;
exit;